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Abstract. An efficient gradient-based method to solve the volume constrained 
topology optimization problems is presented. Each iterate of this algorithm 
is obtained by the projection of a Barzilai-Borwein step onto the feasible set 
consisting of box and one linear constraints (volume constraint). To ensure the 
global convergence, an adaptive nonmonotone line search is performed along 
the direction that is given by the current and projection point. The adaptive 
cyclic reuse of the Barzilai-Borwein step is applied as the initial stepsize. The 
minimum memory requirement, the guaranteed convergence property, and al- 
most only one function and gradient evaluations per iteration make this new 
method very attractive within common alternative methods to solve large-scale 
optimal design problems. Efficiency and feasibility of the presented method 
are supported by numerical experiments. 

Keywords and phrases. Barzilai-Borwein step-size, distributed parameter 
identification, large-scale topology optimization, method of moving asymptotic 
(MMA), nonmonotone line search, volume constraint. 



1. Introduction 

The goal of topology optimization is to find optimal material distribution in 
the given design domain subject to some constraints governed by certain physical 
properties and/or some other practical constraints during the design. In the past 
two decades, advances in the theory of homogenization, optimization, numerical 
analysis as well as newly developed engineering approaches make topology opti- 
mization techniques to become a standard tool of engineering design, in particular 
in the field of structural mechanics. For more detailed literature review, one may 
refer [1, 3] and the references therein. In topology optimization, the design param- 
eters are often material properties (e.g., conductivity or stiffness tensor) and the 
objective is often to minimize an integral functional defined on the spatial domain 
with state variables satisfying a partial differential equation (PDE) corresponding 
to certain physical law. In addition, some bound constraints on the design variables 
and usually a global material resource constraint are also often enforced during the 
design. 
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One typical large class of topology optimization problems have the structure 
that a nonlinear non-convex objective functional is minimized over a feasible region 
defined by a second order elliptic PDE together with bilateral bound and a single 
equality constraints. In this paper, we focus on solving this class of optimal de- 
sign problems which have a wide range of applications in engineering design, such 
as compliance mechanics, fluid dynamics, heat transfer, functionally graded and 
composite materials [1, 3]. 

In the nonlinear optimization literature, there are many well-known optimization 
methods for finding solutions of general optimization problem. Among them, the 
sequential quadratic programming (SQP) [14] method is widely used and generally 
considered to be an efficient method for smooth nonlinear optimization problems 
with constraints. Trust region methods [7] are another class of well-studied meth- 
ods which have strong global and local convergence properties. More recently, 
interior-point methods receives much more attention because of their polynomial 
complexity. However, most of the above mentioned methods require the evaluation 
or a certain type of approximations of the hessian at each iteration. In addition, 
these methods also often require to solve a linear system of equations (usually in- 
definite, dense and is difficult to solve by iterative methods) to a certain level at 
each iteration. Therefore, when the problem size is very large, which often occurs 
after the discretization of the topology optimization problem, obtaining the Hessian 
information as well as solving very large linear system of equations at each iteration 
could be very expensive. Hence, although those second order methods enjoy fast 
local convergence properties, they are generally not efficient to solve very large- 
scale problems, especially when a solution with very high accuracy is not strictly 
required. 

On the other hand, some classes of optimization algorithms are developed within 
engineering community to solve large-scale engineering design problems. The first 
method of this kind is known as CONLIN or convex linearization [12]. More ad- 
vanced version of CONLIN, called method of moving asymptotes (MMA), was in- 
troduced by Svanberg [19] . Within each iteration of this method the optimization 
problem is approximated by a convex separable sub-problem, for which efficient 
solvers are available. The globalization of the method is performed by either line- 
search [21] procedures or conservative sub- iterations [20]. 

Despite the improvements of these methods for general structural optimization 
problems, a key issue to design an efficient optimization method is to exploit the 
specific structure of the desired problem. The optimality criteria (OC) method [3] 
is the first method developed particularly to solve the resource constrained topology 
optimization problems. Recently, OC becomes very popular and is the most widely 
used method in the engineering community for such problems. However, OC is 
not globally convergent and its application is limited to some special problems like 
thermal/structural compliance minimization [?, see:]ch.5]allaire2002soh. Moreover, 
the convergence rate of OC is not very promising. 

By applying the recent techniques developed in the field of nonlinear optimiza- 
tion, in this paper we would like to design an optimization algorithm for solving 
very large-scale topology optimization problems. Each iterate of this algorithm is 
obtained by performing an adaptive nonmonotone line search along the line segment 
connected by the current point and the projection point of a cyclic Barzilai-Borwein 
(CBB) step onto the feasible set. Hence, the method can be called the projected 
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cyclic Barzilai-Borwein (PCBB) method. Because of the special structure of the 
feasible set of the problem, which consists of box constants and a single linear con- 
straint, it is possible to do the projection on the feasible set very efficiently with 
linear time complexity. The main attractive features of the presented algorithm are: 
producing strictly feasible iterations; using almost one objective function and gra- 
dient evaluations per iteration; O(n) memory consumption (6n working memory); 
easily to be implemented; only first order (gradient) information being required. 

2. Barzilai-Borwein methods 

Consider to solve the following finite-dimensional unconstrained optimization 
problem, 

(2.1) min/(x), x£8", 

where / : R™ — >• R is continuously differcntiable, and R™ denotes the Euclidean 
space with dimension n. Suppose that Xo is the starting point, X& is the current 
point, and g^ is the gradient of / at x^, i.e., g k = V/(xfc). Then gradient methods 
generate the next iterative point by 

(2.2) x I+ i = x fe - a k g k , ft = 0,1,..., 

where the stepsize a k is computed by some line search techniques. Two classical 
ways of selecting initial stepsize in the line searches are given by the so called Steep- 
est Descent (SD) and Minimal Gradient (MG) methods, which minimize f(x k —ag k ) 
and ||g(xfc — agfc)|| along the search direction g fc , respectively: 

(2.3) al D = argmin/(x fe - ag fe ), 
and 

(2.4) af G = argmin||g(x fe -ag fe )||. 

where || ■ || denotes the Euclidean norm of a vector. However, it is well-known 
that SD and MG methods can be very slow when the Hessian of / is singular or 
nearly singular at the local minimum. In this case the iterates could approach the 
minimum very slowly in a zigzag fashion [13]. 

The basic idea of Barzilai-Borwein (BB) [2] method is to use the matrix D(a&) = 
^1, where I denotes the identity matrix, to approximate of the Hessian V 2 /(xfc) 
by imposing a quasi-Newton condition on D (of/.): 

(2.5) al B = argmin ||D(o) s k ^ - y fc _i|| 2 , 

where s fe _i = x fc - x ft _ x , y fc _i = g fc - g fe _!, and k > 2. By straightforward 
calculation, BB stepsize obtained from (2.5) is 

(2.6) of* = ^t±. 

Bfc-iyfc-i 

By symmetry, another alternative BB stepsize could be computed by: 

(2.7) af 52 = arg min \\ Bk -i - T> 1 (a) y fc _x || 2 , 

which gives 
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In contrast to the SD or MG methods in which the non-trivial and expensive opti- 
mization problem (2.3) or (2.4) need to be solved to obtain the initial stepsize, the 
BB stepsize is readily available during the iterations by formula (2.6) or (2.8). In 
practice, to keep the stability of the numerical procedure, it is often to project the 
BB stepsize onto a safeguard interval, that is to set 

(2.9) af B = rnm{a max , max{a rara , af B }} 

where a mm ,a max e R and < a min « 1 << a max < oo. 

The following Lemma shows the spectral property of Barzilai-Borwein methods 
and it is due to this property that these methods are usually called spectral gradient 
methods. 

Lemma 2.1. The Barzilai-Borwein stepsize, a BB , is the inverse of the Rayleigh 
quotient, related to vector Sk—i, of the averaged Hessian of the objective function 
between two consecutive iterations k — 1 and k. 

Proof. By the Mean- Value Theorem and straightforward computations, one has 
Sfe - gfc-i 



X fe - X fc _! 



/ V 2 /(ix fc + [1 - t]xfc_i)£ft. 
Jo 



Therefore, 



s t iyfc _ i s^_ 1 ^/ 1 V 2 /(xA ; -i+is fe _i)dtJs fe _i 

s k-l S k-l s k-l S k-l 

which complete the proof. □ 

By Lemma 2.1, one has A mm (af ^ A max , where A. max and A TO j ra arc the 
maximum and minimum eigenvalues of the averaged Hessian matrix V 2 /(x^ + 
tSk) dt respectively. Therefore, a BB I can be considered as an approximation of 
the inverse of the averaged Hessian of the objective function. This shows that the 
Barzilai-Borwein methods could incorporate certain useful second order informa- 
tion with extremely little additional expense of computational and memory cost 
compared with SD or MG methods. Due to its easy implementation, efficiency 
and low storage requirement, BB-type methods have been widely used in many 
applications. Exceptional good performances of BB-type methods have been ob- 
served for solving large-scale problems in particular when only approximate (not 
very accurate) solutions are desired, cf. [11]. 

It has been shown that if the exact steepest descent step (2.3) is reused in a cyclic 
fashion, the convergence speed of gradient methods can be greatly accelerated. 
However, it is often very expensive or impractical to find the exact stepsize along 
the steepest descent searching direction for large dimensional problems, unless the 
objective function is a quadratic function. Hence, for non-quadratic objective func- 
tions, it is unrealistic to apply some cyclic fashion of the steepest descent method. 
On the contrary, the BB stepsize (2.6) or (2.8) can still be easily calculated. Hence, 
analogous to the cyclic steepest descent method, the cyclic BB (CBB) has been 
introduced in [9] for general nonlinear optimization. The numerical results given in 
[9] show the CBB methods have excellent numerical performances compared with 
SD methods and even be competitive to some well-known nonlinear conjugate gra- 
dient methods. Given an integer m ^ 1, which is the cycle length, CBB stepsize 
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can be expressed as 

(2.10) agf£=a*f +1 fori = l,...,m, J = 0,1,..., 

The R-linear convergence of CBB method for a strongly convex quadratic objec- 
tive function has been proved in [8], while the local R- linear convergence for the 
CBB method at a local minimizcr for general nonlinear objective function has been 
established in [9]. 

3. Projected Barzilai-Borwein methods 
Consider the following constrained counterpart of problem (2.1) 

(3.1) min/(x), x 6 V, 

where V C W 1 is a closed non-empty convex set. Because of the constraints in (3.1), 
the iterations generated by (2.2) may lie outside of the feasible set V. Therefore, 

(2.2) need to be modified in order to maintain the feasibility of iterations. Gradient 
projection methods (sec: [18]) keep the feasibility of iterates by frequently project- 
ing trial steps generated from (2.2) onto the feasible set. Considering x as the trial 
step, the projection point y of x onto T>, denoted by 'P-d[x], can be computed by 
solving the following minimization problem 

(3.2) y := V v [x] = argmin -||x - z\\ 2 2 . 

zGX> Z 

Since the feasible set is convex, problem (3.2) always has an unique solution. How- 
ever, in general (3.2) is still a convex constrained quadratic programming problem, 
which could be as difficult as the original problem. And for general convex con- 
strained large-scale problems, this projection step at each iteration could be very 
time consuming and is normally the most expensive part of gradient projection 
methods. Hence, there is little interests on applying the gradient projection meth- 
ods for large-scale problem unless the gradient projection step can be performed 
very efficiently. However, in some special cases when efficient algorithms for cal- 
culating the projection (3.2) exist, for example there are only box or single ball 
constraints, these gradient projection methods will be attractive. Fortunately, as 
we will see in the next section, the projection step can be performed very efficiently 
for the volume constrained topology optimization problems. This makes it pos- 
sible for us to design gradient projection type algorithms for volume constrained 
topology optimization. 

Combining the Barzilai-Borwein stcpsizc rule and the gradient projection method, 
the projected Barzilai-Borwein (PBB) method was first introduced in [5]. In this 
method the iteration updating formula (2.2) is modified to 

(3.3) x x+1 =x fc +/? fc d£, fc = 0,l,..., 

where f3 G K + and the search direction d£ (descent direction) is computed by 
connecting the current point to the projection point of the trial iterate (2.2) based 
on the BB stepsize, that is 

(3.4) d< k *=V v [x k -a k 3B g k }-x k . 

In [5], d£ was called spectral projected gradient. It is not difficult to show that 
d k is a descent direction (see Lemma 3.1). This together with the convexity of the 
feasible set V would imply that for a sufficiently small an iterate of (3.3) will 
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reduce the objective function value while simultaneously preserve the feasibility of 
the iterates. 

Lemma 3.1. for all % 6D and a k B > 0. 

(i) (&(*),<£(*)) <^ K(*)|| 2 . 

k 

(ii) ($k(x) — if and only if x is a stationary point for (3.1). 
Proof, see the Proposition 2.1 in [16]. □ 

4. Globalization by nonmonotone line search 

Using descent directions and simply applying the SD, MG or BB stepsizes in 
gradient projection methods is generally not sufficient to ensure global convergence 
of the iterates starting from an arbitrary initial point. Hence, to deal with general 
nonlinear objective function, a globalization strategy is required to guarantee the 
global convergence of the algorithm. 

Monotonic gradient-based methods usually generate a sequence of iterates for 
which a sufficient decrease in the objective function (or the related merit function) 
is enforced at every iteration. In many cases, the globalization strategy accepts the 
stepsize in the search direction, if it satisfy the well-known Wolfe or Armijo type 
conditions (cf. [?, ]ch. 3]noccdal2006no). This can be accomplished using either of 
monotonic line searches or trust region methods. 

Since the search direction is parallel to the negative direction of the gradient 
(projected gradient) in BB (PBB) methods, the globalization of BB (PBB) meth- 
ods by monotonic function value reduction often reduces them to the classic SD 
(projected SD) method. Hence, these monotonic methods will often distroy all of 
the advantages of BB (PBB) methods in contrast to SD (projected SD) method. 
To maintain the inherit spirits of BB-type methods, it is essential to accept the 
initial BB-type stepsize as frequently as possible while simultaneously ensure the 
global convergence. Hence, some nonmonotone line search techniques need to be 
developed to globalize the BB-type methods. The first nonmonotone line search 
technique was developed in [15] in which the main goal was to accept the full New- 
ton step as much as possible. Combing the type of nonmonotone line search in [15], 
the first globalized version of BB nethod (GBB) was introduced in [17]. Following 
[17], the globalized PBB method (GPBB) was suggested in [5]. In these meth- 
ods the following (weaker) objective function value decrease condition is enforced 
during each iteration 

(4-1) /(x fc+ i)< max f(x k -j) + S gl d k , 

where S G (0, 1), is the search direction and m k is a nonnegative nondecreasing 
integer, bounded by some fixed integer M . More precisely 

mo = and ^ ^ min{mfc_i + 1, M} for k > 0. 

Based on the same motivations, the more efficient and adaptive nonmonotone line 
searches were particularly designed for BB-type methods in [16]. The globalized 
projected cyclic Barzilai-Borwein (PCBB) algorithm based on these new nonmono- 
tone line search techniques can be described as follows: 

Algorithm 4.1. 
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Parameters: 

• e G [0,oo), error tolerance. 

• 5 G (0, 1), descent parameter used in Armijo line search. 

• rj G (0, 1), decay factor for stepsize in Armijo line search. 

• ct m i n ,a ma x 6 (0,oo), safeguarding interval for BB stepsize. 

Initialization: 

• k = 0, Xq = starting guess, and f r _ x = f(xo). 

Main Loop: While \\Vv[x k — g k ] — Xk\\oo > £ 

1. Choose ctk G [a mm ,a ma x}. 

2. Compute dk = Vv[xk - akg k ] - Xk- 

3. Choose f r k such that /(a*) < /J ^ max^-i./** *} /£ < J*™* 
infinitely often, where f™ ax — max{/(3fe_i) : ^ i ^ min(fc, M — 1)}. 

4. Lerj/ fl 6e either f£ ormmj/J,/™ 1 }. 

5. Nonmonotone line search: 

5.1. If f(x k + d k )^f R + 8 g T k d k thenp k = l. 

5.2. _E7.se = r/ 3 , where j > is the smallest integer such that 
fix, ■ ,!*<!,.:■ f l: ■ rfS g^d fc . 

6. 5ei iC/t+i = ccfc + /3fc^fe fc = fc + 1. 
£W Main Loop. 



The variable /£ in Algorithm 4.1 denotes the so called "reference" function value 
in nonmonotone line search. It can be seen that the traditional monotone line search 
simply corresponds to the choice of setting f k = /(xfe) at each iteration. And the 
nonmonotone line search developed in [15] corresponds to the choice of setting 
fk = /™ aa: . In our present study, f k is chosen based on Algorithm 4.2 adapted 
from [16]. Let f k denote /(x^). In the algorithm 4.2, the integer a counts the 
number of consecutive iterations for which j3 k = 1 in Algorithm 4.1 is accepted and 
the Armijo line search in step 5 is skipped. The integer I counts the number of 
iterations since the function value is strictly decreased by an amount A > 0. 

Algorithm 4.2. 



R0. If fc = 0, choose parameters: A > L > 0, 71, 72 > 1, and A > 0; initialize 

a — I — U and J — / — / — J_ 1 — J Q . 

Rl . Update f k as follows: 

Rl.l. Ifl = L, then set I = 0, and 

{fmaxmin ■ f -~> „, 

f k /r™-/r n * 71 ' 

f™ ax otherwise. 

R1.2. Else If a> A, then set 



Jk % J h > Jk and /™"*-/ fc ^ 72, 

fk-i otherwise. 



R1.3. Else set f r k = f^_ v 
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R2. Set f R as follows in step 4 of Algorithm 4.1: 

R2.1. If j = (the first iterate in a CBB cycle) then f R = fT. 
R2.2. Else (3 > 0) f R = min{/£, 

R3. If Bk = 1 in Algorithm 4-1 then a = a + 1, Else (j3k < I) a = 0. 

R4. // f k+1 ^ ff in - A then set f™*™™ = f™ft = f k+1 and I = 0; Else set 
1 = 1 + 1, = fjp in and f™* min = max{fJT axmm , f k+1 }. 

The variable f k nax in Algorithm 4.2 stores the maximum of recent function val- 
ues and /™ m stores the minimum function value within the tolerance A. The 
variable fjj laxmm stores the maximum function value since the last new minimum 
was recorded in f k Tim . 

The condition / (xfc) < in step 3 of Algorithm 4.1 guarantees that the Armijo 
line search in step 5 can be satisfied. Notice that the requirement "/£ < /™ ax infin- 
itely often" in step 3 which is required to ensure the global convergence is a weaker 
condition. Besides Algorithm 4.2 which satisfies this condition, this condition can 
be satisfied by many other strategies. For example, it is possible to set /£ = /™ ax 
at every L iterations. In Algorithm 4.2, f r k = /™ aa: if f( x k-L) - /(x fc ) < A for 
given decrease parameter A > and integer L > 0. 

Now lets give more details about computation of a k in step 1 of Algorithm 4.1. 
This parameter is computed based on the safeguarded CBB scheme as follows. Let 
j as an integer that counts the number of times in which the current BB step has 
been reused and let m as the CBB memory in (2.10), i.e., the maximum number of 
times the BB step will be reused. 

Algorithm 4.3. 

50. If k = choose ao £ [ctmin, OL max \ and a parameter 6 < 1 near 1; set j = 
and flag=l. If k > set flag = 0. 

51. < \dki\ < otk\gki\ for some i (component of vector) then set flag = 1. 

52. // 8k = 1 in Algorithm 4-1 then set j=j+l; Else (0k < i) set flag =1. 

53. If j > m or flag=l or sly k /\\s k \\\\y k \\ ^ 6 then: 
S3.1. If sly k < then 

53. 1.1. If j > 1.5m then set t = min{||a:fc||oo, l}/||d 1 (a;fc)||oo, 
a k+ i = min{a max ,max{a m i„,^ fc }} and j = 0; 
where d 1 ^) = Vv[x k - g k ] - x k 

53. 1.2. Else afc+i = a k 

S3. 2 Else set a k+ i = min{a max , max{a m i„, ce k B }} and j = 0. 

In Algorithm 4.3, the former BB stepsize is reused for the current iterate unless 
one of the following conditions happens in which the new BB stepsize is computed 
(see S3. 2) : (I) the previous BB stepsize is truncated by the projection step, i.e., 
when the trial point using BB stepsize lies outside of the feasible domain and 
the gradient projection was performed (see SI); (II) the previous BB stepsize is 
truncated by the line search step (see S2 where k < 1); (III) the number of 
times the BB stepsize was reused reaches to its bound, i.e., j ^ m (see S3); (IV) 
s fc yfc/ll s fc|| lly/JI i s cl° se to 1 (see [?, ] section 4]dai2006cbb for details about the 
justification for this decision). The condition s^y fc < (see S3.1) is equivalent to 
detection of the negative curvature in the searching direction. Assuming that the 
objective function can be well approximated by a quadratic function in the vicinity 
of the current iterate, a relatively large stepsize should be used in the next iteration 
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(see S3. 1.1) to reduce the function as much as possible once a negative curvature 
is detected. This strategy is similar to the original SPG algorithm (see [?, Jsection 
2]birgin2001ass). 

Now lets briefly review the convergence theory of Algorithm 4.1. 

Theorem 4.4. Let C be the level set defined by 

C = {x g V : /(x) < /(x ) }. 
We assume the following conditions hold: 

Gl. / is bounded from below in C and d max = sup fe ||dfc|| < oo. 

G2. // C is the collection of x 6 T> whose distance to C is at most d max , then 
V/ is Lipschitz continuous on C 
Then either algorithm 4-1 with e = terminates in a finite number of iterations at 
a stationary point, or we have liminf ||d 1 (a;fc)|| 00 = 0. 

k— > oo 

Proof, see the Theorem 2.2 in [16]. □ 

When / is a strongly convex function, (3.1) has a unique minimizer x* and the 
conclusion of the global convergence Theorem 4.4 can be strengthened as follows. 

Theorem 4.5. Suppose f is strongly convex and twice continuously differ entiable 
on T>, and there is a positive integer L with the property that for each k, there exists 
j G [k,k + L) such that fj < /j nax . Then the iterates Xfc of Algorithm J^.l with 
e = converge to the global minimizer x*. 

Proof, see the Corollary 2.3 in [16]. □ 




Figure 1. Projection onto the feasible set for n = 2 related to the 
volume constrained topology optimization problem; the feasible set 
is a potion of line a T x = b which is located inside the box 1 ^ x ^ 
u. 



5. Projection onto the feasible set 

The introduced algorithm in the previous section is general and can be applied 
to any type of convex constrained optimization problems which satisfy the require- 
ments of Theorem 4.4. However, the performance of the method is significantly 
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affected by the efficiency of the projection step. In this section, we explore the spe- 
cial structure of the feasible set in the volume constrained topology optimization 
and introduce a very efficient algorithm to do the projection step. 

After discretization, the feasible set in the volume constrained topology opti- 
mization problems has the following form (see Figure 1) 

(5.1) P={xe K n : a T x = b, Ux^u}, 

where a e f", a, e t +! 6 e I + . l,u e K". < (, ^ ii, < oo. Henceforth, in this 
section we denote by T> the feasible domain defined in (5.1). 

Considering (3.2) together with (5.1), for a given trial vector x, 7 ? x>[x] is the 
unique minimizcr of the following box constrained Lagrangian: 

(5.2) £(z;A) = i||z||2-x T z + i||x||2 + A(a T z-6), z E B, 

where A e 1 is a proper Lagrange multiplier related to the volume constraint and 
the box constrained set B is defined as B = {z E M™ : 1 ^ z ^ u}. For any fixed 
value of A, (5.2) is a convex separable minimization problem with respect to z, 
which has the following explicit solution: 

(5.3) z(A) = max{l, min{x — Aa, u}}, 

where the max and min operators are understood as componentwise. Since solutions 
resulted from (5.3) satisfy the bound constraints, to solve the projection problem 
(3.2) the remaining task then turns out to find the Lagrange multiplier A* such 
that (5.3) satisfies the volume constraint. Hence, an efficient algorithm is needed 
to find the (unique) root A* of the following nonlinear non-smooth one-dimensional 
equation: 

(5.4) fl (A) = a T z(A) -6 = 0. 

Considering (5.3) together with (5.4), it can be seen that the graph of g(X) has 2n 
breakpoints at: 

Aj = (xi — k) j Oi, A" = (xi — Ui)/a,i, i = l,...,n. 

Since k ^ Ui and a,i > for all i, we have A" ^ A'. So, each Zi(X) in (5.3) can be 
expressed in the following form: 

f u h if As^A' 1 , 

(5.5) Zi{X) = l Xi-Xai, if XV^X^Xl 

{ k, if X^X\. 

From (5.5), it is obvious that for each z, Zi{X) is a continuous piecewise linear and 
non-increasing function of A (sec Figure 2). Therefore, by ai > for all i, g(X) is 
a continuous piecewise linear and non-increasing function of A. Then, with some 
straitforward calculations, we can see that the root A* of g(X) is unique and A* E 
[A m m; A maa; ] with X m i n ^ ^ X max: where A m j n = min{A^} and X max = max{A'}. 
So, starting from interval [A m i„, X max ] and using some classical one dimensional root 
finding method, it is possible to find A* up to the machine precision with a priori 
known computational complexity. However, it is possible to use more advanced 
root founding methods to improve the performance of this step. 

In our approach, the Brent's root finding method [6] is employed to solve (5.4). 
The Brent's method does not assume the function differentiability and is enable to 
manage the limited precision of the computed arithmetic very well. In the worst 
condition, the convergence of this method is never slower that of the bisection 
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method. The Brent's method has been proved to be very efficient and robust in 
practice and it is currently accepted as a standard method for one dimensional 
root finding problem (cf. [?, ] chapter 9]press2007nre). The specific implementation 
details of Brent's method are available in the Numerical Recipe (see [?, Jchaptcr 
9]press2007nre).. 
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Figure 2. Plot of Zj(A) as a function of A, cf. equation 5.5. 



6. Numerical solution of topology optimization problems 

To show efficiency of the proposed method for the volume constrained topology 
optimization [4], in this section we solve the standard model problems presented 
in [10]. The description of the problem is as follows: considering two isotropic 
conducting materials, with thermal conductivities k a and hp, < k a < hp, in 
a simply connected design domain ft C M. d (d = 2,3), the goal is to mix these 
materials with a fixed ratio in ft, such that the total temperature gradient in ft is 
minimized under the thermal load / £ L 2 (ft). More specifically, the problem can 
be formulated as 



argmnXy, J(w) = 
subject to: 

-v • oo)w) = 

0(x) = 
k(w) = 

Jn wdx = 



i/ n V0.V0dx, 



/(x) in ft, 

0o{x) on dtt, 

w p kp + (1 - w p )k a in ft, 

R\Cl\, 0<R<1, 0<iy<l, 



where 8 £ i? 1 (ft) is the state variable, p ^ 1 is a penalization factor and w £ L 2 (ft) 
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is the control parameter (topology indicator field). By some standard derivations, 
the first order optimality condition of problem P can be written as the following: 



(OC) 



-V ■ (fc(iu)V0) = 
0(x) = 
-V • (k(w)Vrj) = 

7?(x) = 

k(w) = 

MG(x)) = 
G(x) 



/(x) 
0o (x) 
-V ■ (V0) 


+ (f - w p )k a 





= —pw p 1 



in O, 

on 90, 

in Q, 

on 90, 



m 
in 



^ - k a )V9 ■ Vrj in 0, 



where 77 G i?Q (O) is the adjoint state, G is the L 2 gradient of the objective functional 
with respect to w and Vx>{u) denotes the L 2 projection of function u onto the 
admissible set P, 



V = {weL 2 {tt)\ f w{x)dx = R\Cl\, 0<R<1, O^w^l}. 
Jn 

By discretization of O into an n control volumes, we have the finite dimensional 
counterpart of problem (P) . We also assume that the state variable and the design 
parameter are defined at the center of each control volume. Under these assump- 
tions, the admissible design domain T> forms a simplex in R n which is identical to 
the continuous knapsack constraints in (5.1). At each iterate of the control pa- 
rameter w, we solve the associated state PDE. Then the discretized optimization 
problem would have the general format of problem (3.1), which has a nonlinear 
objective function with convex continuous knapsack constraints. 

In our numerical experiments, the problem (P) was solved in two and three 
dimensions for O = [0,1] 2 and = [0, l] 3 . The spatial domain is divided into 
127 2 and 31 3 control volumes in two and three dimensions respectively In all of 
experiments, we set k a — 1, /(x) = 1 and R = 0.4. And in these experiments two 
conductivity ratios 2 and 100 were tested, which is equivalent to setting hp = 2, 100 
respectively The penalization factor p is taken to be 1 and 10 for conductivity 
ratio 2 and 100 respectively. The governing PDE are solved by cell centered finite 
volume method using central difference scheme and the related system of linear 
equations are solved by a preconditioned conjugate gradient method with relative 
convergence threshold 10 -20 . The optimization is performed for 15 iterations in 
these numerical experiments. Notice that using finite volume method we do not 
observe any topological instability phenomena (the checkerboard pattern), which 
often occurs by using finite element method for topology optimization problems. 

The input parameters related to optimization algorithms used in this study are 
as follows: 5 = 10~ 4 , 17 = 0.5, a mm = 10~ 30 , a max = 10 30 , A = 40, L = 10, 
M = 20, m = 4, 71 = 72 = 2, 9 = 0.975. Moreover, the initial value ao for spectral 
step-size was taken to be l/UPpfxfe — g fc ] — x^H^. Another alternative choice for 
this parameter could be a = l/H^ffx*; — gfc] — X/blh- However, for the testing- 
problems in our numerical experiments the former choice was considerably more 
efficient. 
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To evaluate the efficiency of the presented method, we compare our results with 
the results obtained by the method of moving asymptotic (MMA) [19]. The imple- 
mentation of MMA available in SCPIP code [22] 1 is used in our experiments. All 
default parameters are used in SCPIP code, except the parameter for the constraint 
violation, for which the threshold 10~ 7 is used in this study. That SCPIP code used 
in MMA algorithm has two globalization strategics. The first one is identical to 
that of the original MMA by Svanberg [20], while the second one is globalization 
by monotone line search method (see: [21]). The first strategy is employed in our 
experiments, since our results with second strategy was significantly worse in terms 
of the computational cost. Note that in our procedure the PDE constraint is solved 
upto the accuracy of the finite volume method we have applied for solving this 
PDE, and the constraints (5.1) in PCBB algorithm are satisfied upto the machine 
precision. 

The results of our numerical experiments including the variation of the objective 
functional during the optimization process and the final resulted topology (w-ficld) 
are shown in figures 3, 4, 5 and 6. The plots in figures 3, 4 show the success 
of the presented method to solve these topology optimization problems. Roughly 
speaking, both methods behave very competitively in terms of computational cost 
and final results. The main differences in the results are related to the conductivity 
ratio 100, in which the PCBB performs superior and its final objective function 
values are considerably lower than that of MMA. The sign of the differences is 
clear in the figures 5 and 6 related to the final topologies. But for the conductivity 
ratio equal to 2, it seems MMA behaves slightly better than PCBB. However, the 
differences in terms of the objective function values as well as the final topologies 
are almost negligible . 

In all our numerical experiments, the final total number of function and gradi- 
ent evaluations was 15 which is equal to the number of optimization cycles. This 
result shows that both methods used only one function and gradient evaluations 
per iteration in practice. We believe that this is a key property for the success of 
MMA and makes the method well accepted in the engineering design community. 
In fact, MMA behaves very conservatively and uses small steps to proceed toward 
a local minimum. More clearly, it does not use (expensive) line search globalization 
strategy (unlike alternative methods), but uses reasonably small steps such that 
hopefully the merit function decreases sufficiently during each iteration. Of course, 
whenever the merit function increases (or sufficient decrease in merit function value 
violates), which rarely happens in practice, it performs sub-cycles to ensure the de- 
sired monotonic behavior. On the other hand, the nonmonotone PCBB methods 
enjoys such property using an alternative strategy. In practice, by applying the 
nonmonotonic line search, PCBB often uses only one function evaluation per opti- 
mization cycle. As our results clearly show, this property makes the nonmonotone 
PCBB method as a very competitive alternative to MMA for these class of prob- 
lems. 

It is important to note that in all of our experiments presented here, the pro- 
jection step is used actively in PCBB method. Therefore, the cyclic reuse of step 
size never employed in our testing problems (cf. algorithm 4.3). Moreover, for the 
conductivity ratio 100, PCBB method explored directions of negative curvature 



The SCPIP code (in Fortran) is freely available through personal request from its original 
author (Christian Zillober: christian. zillobcr(at)uni-wucrzburg.dc). 
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after a few iterations, and so used large step-sizes which considerably accelerate 
the convergence. This property plays an important role for the superior results of 
PCBB compared with MMA in these cases (cf. table 1). 

Besides the presented numerical results here, we have also applied the presented 
method successfully to many families of topology optimization problems with very 
satisfied results. But we omit the details of these experiments here due to the 
limited space for this paper. 




2D mailt PCBB with conductivity rulio = 100 2D result: MMA with conductivity ratio = 100 



Figure 3. Variations of scaled objective function values during 
optimization cycles for 2D examples. 

Table 1. Variations of step size during the fist 10 optimiza- 
tion cycles in PCBB algorithm for 2D examples (the decimals are 
dropped in the table). a(2) and a(100) denote values of oik for the 
conductivity ratio 2 and 100 respectively. 



iter 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


a(2) 


37 


48 


124 


133 


52 


44 


54 


147 


305 


286 


a(100) 


53 


10 3u 


10 3u 


243 


164 


100 


104 


10 au 


10 3U 


1093 
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Figure 4. Variations of scaled objective function values during 
optimization cycles for 3D examples. 




Figure 5. Final distribution of material field, w, for 2D examples 
related to PCBB (top) and MMA (bottom) for conductivity ratio 
2 (left) and 100 (right). 




Figure 6. Final distribution of material field, w, for 3D examples, 
related to PCBB (top) and MMA (bottom) conductivity ratio 2 
(left) and 100 (right). 
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7. Closing remarks 

Recent studies on the spectral projected gradient methods show that these class 
of methods are very promising for solving the large-scale convex constrained opti- 
mization problems when the projection on the feasible set can be performed effi- 
ciently. In this paper, we presented a particular spectral projected gradient method 
called PCBB (Projected Cyclic Barzilai-Borwein) for solving volume constrained 
topology optimization problems. This method applies the cyclic Barzilai-Borwein 
stepsizes and uses the most recent adaptive nonmonotone line search techniques, 
which greatly improves the efficiency of the method as well as ensures its global 
convergence. By exploring the structure of the admissible set of the volume con- 
strained topology optimization problem, the projection step can be performed very 
efficiently. In addition to high efficiency, our presented method also enjoys the fol- 
lowing features: easy implementation, minimum memory requirement, no need for 
second order (Hessian) information, not sensitive to data noise, feasibility being 
strictly maintained during optimization process. All these features are essential for 
a successful numerical method to solve large-scale topology optimization problems. 

Our numerical results indicate our presented methods are very promising and 
well-suited for the class of topology optimization problems considered in this paper. 
Comparing our results with those of MMA, a well accepted method in shape and 
topology optimization community, it seems the presented methods could be a very 
competitive alternative choice for solving the class of topology design problems. 
Finally, our results suggest that including spectral step size and a nonmonotone 
globalization strategy in MMA, the MMA algorithm could be further significantly 
improved. This would be a topic for our continuing research. 
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